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During  the  period  covered  by  the  Grants  attention  has  been  focused  on 
three  areas,  all  of  them  of  importance  in  the  successful  application  of 
implicit  algorithms  to  Computational  Fluid  Dynamics  (CFD) : 

1,  The  role  of  boundary  concitions  for  implicit  hyperbolic  schemes. 

2.  The  stability  of  hyperbolic  Approximate  Factorization  schemes 

in  three  space  dimensions,  /'''  '  ■ 

\ 

3i  The  rate  of  convergence  to  steady  state  of  ADI  methods. 

This  report  delineates  the  progress  in  each  of  the  above  enumerated 
areas.  The  details  of  the  research  will  be  found  in  reports  and  papers  as 
referenced  below  for  each  of  the  tasks. 

1.  Role  of  Boundary  Conditions 

It  is  well  known  that  for  problems  governed  by  hyperbolic  partial 
differential  equations  one  cannot  specify  boundary  conditions  along  all  the 
Qomain.  In  fact,  every  point  on  the  boundary  which  receives  infcrmation, 
via  the  characteristics,  from  the  initial  data  must  be  excluded  from 
boundary  val.ue  specification.  On  the  other  hand,  in  the  numerical  solution 
of  the  approximating  difference  equations  one  must  often  use  those  "excluded" 
boundary  points  in  order  to  advance  the  integration  temporally.  (This  may 
be  avoided  by  using  upwind  differencing;  one  then,  hou’ever,  may  pay  a 
price — e.g.  lower  spatial  accuracy.)  The  theory  of  how  to  treat  numerically 
these  "forbidden"  hyperbolic  boundaries  is  due  to  Kriess  et  al  [1,2, 2, 4, 5] 
and  Osher  [6,7],  The  theory  is  not  easy  to  apply  to  "systems  of  equations 
or  in  more  than  one  space  dimension.  In  the  multidimensional  case  one  must 
assume  periodic  solutions  in  the  directions  not  normal  to  the  boundary  being 
investigated.  Application  of  the  G-K-S  theory  to  one-dimensional  Beam.- 
Warming  type  algorithms  was  carried  out  by  Gustafsson  and  Oliger  [8] ,  and 
Yee,  Beam  and  Warming  [9]  for  extrapolation-type  bour;dary  concitions. 

In  the  present  research  it  was  thought  interesting  to  find  out  whether 
the  conclusions  drawn  from  the  one-dimensional  analysis  carry  over  tc  the 
2-D  case.  Some  previous  work  [10]  with  explicit  schemes  indicated  that  this 
is  not  always  the  case.  The  details  of  the  present  investigation  were  reported 
at  the  NASA-Ames  workshcp  on  boundary'  conditions  (Oct.  1981) ,  have  come  out 
as  an  ICASE  report  and  were  published  in  the  JCP  [11].  The  major  conclusions 
of  this  study  may  be  summarized  as  follows:  extrapolation-type  boundary 


2 


conditions  which  are  stable  in  one  dimension  will  rerruair.  i:.  nwc  dimensions 
if  the  indicated  extrapolation  is  carried  out  in  a  direct icr.  nc ritial  tc  the 
plane  of  the  boundary.  If  the  extrapolation  is  dene  in  a  ncn-ncrmal 
direction  (skewed  to  the  boundary  plane) ,  then  the  cutceme  depends  on 
the  nature  of  the  basic  internal-points  algorithm.  If  this  internal  scheme 
is  strictly  dissipative  (e.g.  iMacCormack '  s  algorithr.)  ,  then,  the  1-D  results 
carry  over  even  if  the  extrapolation  in  2-D  be  performed  in  a  direction 
skewed  to  the  boundary  plane.  If,  on  the  ether  hand,  the  internal  scheme 
is  not  strictly  dissipative  (e.g.  backward  Euler) ,  then  the  1-E  boundary 
treatment  remains  stable  only  if  the  extrapolation  is  normal  tc  the  bc-ondary 
plane;  otherwise  it  is  unstable  (not  ever,  conditionally  Fzahle)  . 

2.  Stability  of  3-D  Hyperbolic  Approximate  Faotcrioat in-,  schemes 

It  has  been  pointed  out  by  Dwoyer  and  Thames  [12]  that  the  Beam-Karmiinc 
algorithm  in  3-D  is  unstable  for  the  linear  wave  equation.  It  was  suspected 
but  not  demonstrated  that  this  is  also  true  of  the  system  of  linearized 
equations  resulting  from  Euler's  equations  of  gas-dynamics. 

As  part  of  the  present  research  effort  it  was  undertaken  first  to 
demonstrate  the  instability  of  the  Approximate  Factcritaticr; ,  backward 
Euler  algorithm  when  applied  to  the  equations  of  gas-dyr.anios  (Euler 
equations)  .  The  second  task  was  to  try  and  construe^  an  algorithm',  which 
will  retain  the  advantages  of  the  Beam-Warm.ing  algorithm  (approximate 
f acocrization  of  the  implicit  terms,  "delta  form,"  of  the  unknown  vector)  , 
or  of  the  LODQ  algorithm  (see  [12]),  and  yet  become  stable  due  tc  changes 
introduced  into  the  explicit  terms  only.  Both  tasks  were  completed 
successfully  and  are  reported  in  [13]. 

The  new  stabilized  algorithm  has  the  fcllowir.a  features:  three  LCD 
(locally  one  dimensional)  fractional  steps  are  folicwed  by  pliczt 

step.  The  total  algorithm  is  stable  in  3  dimensicns,  is  spat,  .ly  second 
order  accurate  for  a  large  range  of  the  Courant  number  when  it  converces 
to  steady  state.  The  intermediate  steps  leading  tc  the  steady  state  are 
not  time  consistent  and  hence  do  not  represent  a  truly  evolving  solutio.n. 

This  new  algorithm  when  applied  to  the  2-dimer.sional  case  (which  is  stable 
for  the  standard  AF  methods)  accelerates  the  convergence  tc  steady  state 
by  factors  between  2  to  3  depending  on  the  mesh  size.  The  numerical 
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experiments  carried  out  in  connection  with  this  work  brcucnt  out  the  little 
known  (or  rather  little  acknowledged)  fact  that  the  rate  of  convergence 
of  hyperbolic  AF  algorithms  to  steady  state  is  very  sensitive  to  the  time 
step.  In  fact,  for  the  linear  wave  equation,  the  optimal,  Courant  number 
is  of  order  unity.  Similar  results  for  the  full  Euler  equations  were  also 
reported  by  Thompkins  and  Bush  [14].  These  empirical  observations  led  us 
to  concentrate  on  the  third  topic. 

3.  Converoence  Rate  to  Steady  State  of  ADI  Method 

ADI  methods  for  elliptic  partial  differential  equaticns  were  proposed 
already  in  the  1950' s.  Particularly  v.ell  knov/n  were  the  algorithms  du"^  tc 
Peaceman  and  Rachford  [15]  and  Douglass  and  Gunn  [16].  In  recent  years  ADI 
methods  were  advocated  to  solve  parabolic  and  hyperbolic  partial  differentia 
equations,  see  Beam  and  V.’arrr.ing  [!"’],  and  Briley  and  McDonald  [IS].  The 
motivation  for  this  extension  of  ADI  schemes  was  to  combine  the  convenience 
of  one-dimensional  easily  invertible  operators  and  the  unconditional 
stability  of  implicit  methods.  As  pointed  cut  above,  the  large  expected 
gains  over  direct  explicit  solvers  have  not  been  realized  fully.  This  is 
so  because  the  rate  of  convergence  to  steady  state  is  very  sensitive  tc 
the  Courant  number  and  it  decreases  rapidly,  i.e.  the  iteration  count  grows 
rapidly,  when  the  calculation  is  carried  out  away  from  an  optional  time 
step.  At  the  optimal  Courant  number  the  convergence  rate  is  comparable  to 
those  of  explicit  meth.ods. 

The  work  carried  out  under  tne  present  Grants  concentrated  so  far  cn 
the  parabolic  case  w'hich  is  easier  tJ  analyze.  The  results  obtained  so 
far  may  be  summarized  as  follows  (see  [19])  : 

3.1  The  convergence  to  steady  state  of  parabolic  ADI  solvers,  such 
as  the  Dcuglass-Gunn  or  the  Peacemar.-Rachford  algorithms,  is  analyzed  in 
terms  of  the  norm  of  the  residual.  This  approach,  which  assumes  the 
presence  of  many  frequencies  and  averages  over  their  spectr^um*,  turns  out  to 
be  sudeessful  in  predicting  how  the  number  of  iterations  needed  to  converge 
to  steady  state  depends  on  the  Courant  number. 

3.2  A  new  corrected  ADI  algorithm  has  been  devised  which  has  the 
following  properties: 

a)  Its  construction  necessitates  only  the  addition  of  the  same 
expliri*’  teT-m  to  til  existing  Approximate  Factorization  codes. 
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It  is  robust  in  the  sense  that  it  need  not  be  fine-tuned 
for  different  mesh  sizes,  different  grid  stretchings,  mixed  Dirichlet- 
Neumann  boundary  conditions,  etc. 

c)  The  rate  of  convergence  to  steady  state  is  substantially 
improved  and  is  insensitive  to  the  Courant  number  over  a  large  range. 

d)  Its  method  of  derivation  is  easily  extended  to  the  three- 
dimensional  case. 
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